PCA plots

All code is hidden by default, to show code click on code or select ‘show all code’ on the top right of this document.

Setup

reqpkg <- c("tidyverse","feather", "DESeq2", "dplyr", "magrittr", "genefilter", "AnnotationHub", "org.Mm.eg.db", "GO.db", "vsn", "pheatmap", "biomaRt", "curl", "RColorBrewer", "viridis", "ggplot2", "ggExtra", "ggpubr", "fgsea", "DT", "ggfortify", "lfda", "rgl", "pca3d")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "tidyverse"
## [1] '1.3.1'
## [1] "feather"
## [1] '0.3.5'
## [1] "DESeq2"
## [1] '1.24.0'
## [1] "dplyr"
## [1] '1.0.6'
## [1] "magrittr"
## [1] '2.0.1'
## [1] "genefilter"
## [1] '1.66.0'
## [1] "AnnotationHub"
## [1] '2.16.1'
## [1] "org.Mm.eg.db"
## [1] '3.8.2'
## [1] "GO.db"
## [1] '3.8.2'
## [1] "vsn"
## [1] '3.52.0'
## [1] "pheatmap"
## [1] '1.0.12'
## [1] "biomaRt"
## [1] '2.40.5'
## [1] "curl"
## [1] '4.3'
## [1] "RColorBrewer"
## [1] '1.1.2'
## [1] "viridis"
## [1] '0.5.1'
## [1] "ggplot2"
## [1] '3.3.3'
## [1] "ggExtra"
## [1] '0.9'
## [1] "ggpubr"
## [1] '0.4.0'
## [1] "fgsea"
## [1] '1.10.1'
## [1] "DT"
## [1] '0.17'
## [1] "ggfortify"
## [1] '0.4.11'
## [1] "lfda"
## [1] '1.1.3'
## [1] "rgl"
## [1] '0.104.16'
## [1] "pca3d"
## [1] '0.10.2'

options(rgl.useNULL = TRUE)
# knitr::knit_hooks$set(webgl = hook_webgl)
theme_set(theme_pubr())
select <- dplyr::select
splits <- list("WTSPF"=1:18, "KOSPF"=19:36, "WTGF"=37:54, "KOGF"=55:72)
group.colors <- c("SPF:WT" = "#103EFB", "SPF:KO" = "#FC2A1C", "GF:WT" ="#128D15", "GF:KO" = "#FD9226")
source("ggbiplot.R")
A = windowsFont("Helvetica LT Pro Light")

ensembl <- useMart("ensembl")
ensemblMouse <- useDataset("mmusculus_gene_ensembl",mart=ensembl)
mouseProteinCodingGenes <- getBM(attributes=c("ensembl_gene_id","external_gene_name","description"), filters='biotype', values=c('protein_coding'), mart=ensemblMouse)

All data

Make DESeq object.

filePath = '../../../'
df <- read.csv(paste0(filePath, 'KF_RNASeq_counts_filtered.txt')) %>%
  dplyr::select(-X) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>% 
  set_rownames(.$Geneid) %>% 
  dplyr::select(-Geneid)

gt_list <- rep(c(rep("WT",18), rep("KO",18)),2)
cond_list <- c(rep("SPF",36),rep("GF",36))
time_list <- rep(c(rep("ZT02",3), rep("ZT06",3), rep("ZT10",3), rep("ZT14",3), rep("ZT18",3), rep("ZT22",3)), 4)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list, Time=time_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")

dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)
dds

## class: DESeqDataSet 
## dim: 21847 72 
## metadata(1): version
## assays(1): counts
## rownames(21847): ENSMUSG00000051951 ENSMUSG00000025900 ...
##   ENSMUSG00000095041 ENSMUSG00000095742
## rowData names(0):
## colnames(72): KF_01 KF_02 ... KF_71 KF_72
## colData names(3): Genotype Condition Time

Normalize data.

vsd <- vst(dds, blind = FALSE)
head(assay(vsd[,1:5]), 3)

##                       KF_01    KF_02    KF_03    KF_04    KF_05
## ENSMUSG00000051951 4.406910 4.406910 4.406910 4.406910 4.406910
## ENSMUSG00000025900 4.406910 4.406910 4.406910 4.406910 4.406910
## ENSMUSG00000025902 5.946324 5.921394 5.992112 5.503907 5.774604

colData(vsd)

## DataFrame with 72 rows and 4 columns
##       Genotype Condition     Time        sizeFactor
##       <factor>  <factor> <factor>         <numeric>
## KF_01       WT       SPF     ZT02  1.09295313026657
## KF_02       WT       SPF     ZT02  1.13261609555455
## KF_03       WT       SPF     ZT02  1.02495021156392
## KF_04       WT       SPF     ZT06  1.08800022292013
## KF_05       WT       SPF     ZT06  1.11995649245623
## ...        ...       ...      ...               ...
## KF_68       KO        GF     ZT18 0.893688439524359
## KF_69       KO        GF     ZT18  1.02889018870847
## KF_70       KO        GF     ZT22 0.997131890140167
## KF_71       KO        GF     ZT22  0.89087483855546
## KF_72       KO        GF     ZT22 0.970363030607483

Generate PCA

# calculate the variance for each gene
rv <- rowVars(assay(vsd))

# select the ntop genes by variance
select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]

# perform a PCA on the data in assay(x) for the selected genes
pca <- prcomp(t(assay(vsd)[select,]))

intgroup <- c("Condition","Genotype")
intgroup.df <- as.data.frame(colData(vsd)[, intgroup, drop=FALSE])

group <- if (length(intgroup) > 1) {
  factor(apply( intgroup.df, 1, paste, collapse=":"))
} else {
  colData(object)[[intgroup]]
}
percentVar <- pca$sdev^2 / sum(pca$sdev^2)

# assembly the data for the plot
d <- data.frame(PC1=pca$x[,1], PC2=pca$x[,2], PC3=pca$x[,3], group=group, intgroup.df, name=colnames(vsd)) %>%
  mutate(colors=c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))

3D plots

plot3d(d[,1:3], pch=16, size=10, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()

plot3d(d[,1:3], pch=16, size=10, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), cex=0.7, pos=3, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()

plot3d(d[,1:3], size=0.1, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), col=c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()

PC1 vs PC2

p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
p1

Variations of the plot.

ggMarginal(p1 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)

ggMarginal(p1 + theme(legend.position = "none"), type="density", groupColour = TRUE)

I also used ggbiplot to plot the PCA. It doesn’t recalculate the PCA, it just plots it differently.

p <- ggbiplot(pca, groups = group, var.axes = F, ellipse = T, circle = T, obs.scale = 0.5) + ggtitle("PCA plot via ggbiplot") + scale_color_manual(values = group.colors)
p$layers[[1]] <- NULL
p + geom_point(aes(color = groups), size=3)

PC1 vs PC3

p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
p2

ggMarginal(p2 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)

ggMarginal(p2 + theme(legend.position = "none"), type="density", groupColour = TRUE)

PC2 vs PC3

p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
p3

ggMarginal(p3 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)

ggMarginal(p3 + theme(legend.position = "none"), type="density", groupColour = TRUE)

PC1 vs PC2 vs PC3

s = 2
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors)
ggarrange(p1,p2,p3, nrow=1, common.legend = TRUE)

Time-series PCA

Colors go from lighter to darker.

#DESeq from data set, using data frame for design
dds2 <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                       colData = condition_list,
                                       design = ~ Genotype + Condition + Time + Genotype:Condition+Time)
dds2

## class: DESeqDataSet 
## dim: 21847 72 
## metadata(1): version
## assays(1): counts
## rownames(21847): ENSMUSG00000051951 ENSMUSG00000025900 ...
##   ENSMUSG00000095041 ENSMUSG00000095742
## rowData names(0):
## colnames(72): KF_01 KF_02 ... KF_71 KF_72
## colData names(3): Genotype Condition Time

vsd2 <- vst(dds2, blind = FALSE)
head(assay(vsd2[,1:5]), 3)

##                       KF_01    KF_02    KF_03    KF_04    KF_05
## ENSMUSG00000051951 4.680675 4.680675 4.680675 4.680675 4.680675
## ENSMUSG00000025900 4.680675 4.680675 4.680675 4.680675 4.680675
## ENSMUSG00000025902 6.091255 6.068107 6.133806 5.682343 5.932073

colData(vsd2)

## DataFrame with 72 rows and 4 columns
##       Genotype Condition     Time        sizeFactor
##       <factor>  <factor> <factor>         <numeric>
## KF_01       WT       SPF     ZT02  1.09295313026657
## KF_02       WT       SPF     ZT02  1.13261609555455
## KF_03       WT       SPF     ZT02  1.02495021156392
## KF_04       WT       SPF     ZT06  1.08800022292013
## KF_05       WT       SPF     ZT06  1.11995649245623
## ...        ...       ...      ...               ...
## KF_68       KO        GF     ZT18 0.893688439524359
## KF_69       KO        GF     ZT18  1.02889018870847
## KF_70       KO        GF     ZT22 0.997131890140167
## KF_71       KO        GF     ZT22  0.89087483855546
## KF_72       KO        GF     ZT22 0.970363030607483

rv <- rowVars(assay(vsd2))
select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]
pca2 <- prcomp(t(assay(vsd2)[select,]))
intgroup <- c("Condition","Genotype","Time")
intgroup.df <- as.data.frame(colData(vsd2)[, intgroup, drop=FALSE])
group <- if (length(intgroup) > 1) {
  factor(apply( intgroup.df, 1, paste, collapse=":"))
} else {
  colData(object)[[intgroup]]
}
percentVar <- pca2$sdev^2 / sum( pca2$sdev^2 )
d <- data.frame(PC1=pca2$x[,1], PC2=pca2$x[,2], PC3=pca2$x[,3], group=group, intgroup.df, name=colnames(vsd2))

Color scales

group.colors.t <- c("SPF:WT:ZT02" = "#728EFE", "SPF:WT:ZT06" = "#4A6EFC", "SPF:WT:ZT10" = "#224DFC", "SPF:WT:ZT14" = "#0433F1", "SPF:WT:ZT18" = "#032BC9", "SPF:WT:ZT22" = "#0222A1",
                    "SPF:KO:ZT02" = "#FDA19B", "SPF:KO:ZT06" = "#FD7C72", "SPF:KO:ZT10" = "#FC564A", "SPF:KO:ZT14" = "#FC2A1C", "SPF:KO:ZT18" = "#DD1203", "SPF:KO:ZT22" = "#B50F03",
                    "GF:WT:ZT02" = "#6FEC71", "GF:WT:ZT06" = "#26E329", "GF:WT:ZT10" = "#17B51A", "GF:WT:ZT14" = "#128D15", "GF:WT:ZT18" = "#0E6C10", "GF:WT:ZT22" = "#0A480B",
                    "GF:KO:ZT02" = "#FDAD5D", "GF:KO:ZT06" = "#FD9226", "GF:KO:ZT10" = "#FD850D", "GF:KO:ZT14" = "#DE7002", "GF:KO:ZT18" = "#B65C02", "GF:KO:ZT22" = "#8D4701")
group.colors.t2 <- unlist(lapply(group.colors.t, function(x) rep(x, 3)))
group.colors.t

## SPF:WT:ZT02 SPF:WT:ZT06 SPF:WT:ZT10 SPF:WT:ZT14 SPF:WT:ZT18 SPF:WT:ZT22 
##   "#728EFE"   "#4A6EFC"   "#224DFC"   "#0433F1"   "#032BC9"   "#0222A1" 
## SPF:KO:ZT02 SPF:KO:ZT06 SPF:KO:ZT10 SPF:KO:ZT14 SPF:KO:ZT18 SPF:KO:ZT22 
##   "#FDA19B"   "#FD7C72"   "#FC564A"   "#FC2A1C"   "#DD1203"   "#B50F03" 
##  GF:WT:ZT02  GF:WT:ZT06  GF:WT:ZT10  GF:WT:ZT14  GF:WT:ZT18  GF:WT:ZT22 
##   "#6FEC71"   "#26E329"   "#17B51A"   "#128D15"   "#0E6C10"   "#0A480B" 
##  GF:KO:ZT02  GF:KO:ZT06  GF:KO:ZT10  GF:KO:ZT14  GF:KO:ZT18  GF:KO:ZT22 
##   "#FDAD5D"   "#FD9226"   "#FD850D"   "#DE7002"   "#B65C02"   "#8D4701"

show_col(group.colors.t)

3D plots

plot3d(d[,1:3], pch=16, size=10, col = group.colors.t2)
rglwidget()

plot3d(d[,1:3], pch=16, size=10, col = group.colors.t2)
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), cex=0.7, pos=3, col = group.colors.t2)
rglwidget()

plot3d(d[,1:3], size=0.1, col = group.colors.t2)
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), col=group.colors.t2)
rglwidget()

PC1 vs PC2

p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none") +
  ggtitle("PCA with time gradient")
p1

Variations of the plot.

ggMarginal(p1, type="boxplot", groupFill = TRUE)

ggMarginal(p1, type="density", groupColour = TRUE)

I also used ggbiplot to plot the PCA. It doesn’t recalculate the PCA, it just plots it differently.

p <- ggbiplot(pca, groups = group, var.axes = F, ellipse = T, circle = T, obs.scale = 0.5) + ggtitle("PCA plot via ggbiplot") + scale_color_manual(values = group.colors.t)
p$layers[[1]] <- NULL
p + geom_point(aes(color = groups), size=3)

PC1 vs PC3

p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none")
p2

ggMarginal(p2, type="boxplot", groupFill = TRUE)

ggMarginal(p2, type="density", groupColour = TRUE)

PC2 vs PC3

p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=4) + 
  xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none")
p3

ggMarginal(p3, type="boxplot", groupFill = TRUE)

ggMarginal(p3, type="density", groupColour = TRUE)

PC1 vs PC2 vs PC3

s = 2
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none")
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none")
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=s) + 
  xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
  ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
  coord_fixed() +
  scale_color_manual(values = group.colors.t) +
  theme(legend.position = "none")
ggarrange(p1,p2,p3, nrow=1)

Individual timepoints

filePath = '../../../data/R/by_time/'

df.PCA <- function(data, intgroup, colors) {
  # calculate the variance for each gene
  rv <- rowVars(assay(data))
  
  # select the ntop genes by variance
  select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]
  
  # perform a PCA on the data in assay(x) for the selected genes
  pca2 <- prcomp(t(assay(data)[select,]))
  
  intgroup.df <- as.data.frame(colData(data)[, intgroup, drop=FALSE])
  
  group <- if (length(intgroup) > 1) {
    factor(apply( intgroup.df, 1, paste, collapse=":"))
  } else {
    colData(data)[[intgroup]]
  }
  
  percentVar <- pca2$sdev^2 / sum( pca2$sdev^2 )
  
  # assembly the data for the plot
  d2 <- data.frame(PC1=pca2$x[,1], PC2=pca2$x[,2], group=group, intgroup.df, name=colnames(data))
  
  p <- ggplot(data=d2, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=3) + 
    xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
    ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
    coord_fixed() +
    scale_color_manual(values = colors)
  return(p)
}

df <- read_feather(paste0(filePath, 'time_1_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

gt_list <- rep(c(rep("WT",3), rep("KO",3)),2)
cond_list <- c(rep("SPF",6),rep("GF",6))
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p1 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT2")

df <- read_feather(paste0(filePath, 'time_2_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p2 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT6")

df <- read_feather(paste0(filePath, 'time_3_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p3 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT10")

df <- read_feather(paste0(filePath, 'time_4_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p4 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT14")

df <- read_feather(paste0(filePath, 'time_5_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p5 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT18")

df <- read_feather(paste0(filePath, 'time_6_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  as.data.frame() %>%
  set_rownames(.$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)

vsd <- vst(dds, blind = FALSE)

p6 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT22")
ggarrange(p1,p2,p3,p4,p5,p6, common.legend = T)

Session info

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.1.1                plyr_1.8.6                 
##  [3] pca3d_0.10.2                rgl_0.104.16               
##  [5] lfda_1.1.3                  ggfortify_0.4.11           
##  [7] DT_0.17                     fgsea_1.10.1               
##  [9] Rcpp_1.0.6                  ggpubr_0.4.0               
## [11] ggExtra_0.9                 viridis_0.5.1              
## [13] viridisLite_0.3.0           RColorBrewer_1.1-2         
## [15] curl_4.3                    biomaRt_2.40.5             
## [17] pheatmap_1.0.12             vsn_3.52.0                 
## [19] GO.db_3.8.2                 org.Mm.eg.db_3.8.2         
## [21] AnnotationDbi_1.46.1        AnnotationHub_2.16.1       
## [23] BiocFileCache_1.8.0         dbplyr_2.1.1               
## [25] genefilter_1.66.0           magrittr_2.0.1             
## [27] DESeq2_1.24.0               SummarizedExperiment_1.14.1
## [29] DelayedArray_0.10.0         BiocParallel_1.18.1        
## [31] matrixStats_0.58.0          Biobase_2.44.0             
## [33] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [35] IRanges_2.18.3              S4Vectors_0.22.1           
## [37] BiocGenerics_0.30.0         feather_0.3.5              
## [39] forcats_0.5.1               stringr_1.4.0              
## [41] dplyr_1.0.6                 purrr_0.3.4                
## [43] readr_1.4.0                 tidyr_1.1.3                
## [45] tibble_3.0.6                ggplot2_3.3.3              
## [47] tidyverse_1.3.1            
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.1.4                    tidyselect_1.1.0             
##   [3] RSQLite_2.2.3                 htmlwidgets_1.5.3            
##   [5] munsell_0.5.0                 preprocessCore_1.46.0        
##   [7] miniUI_0.1.1.1                withr_2.4.1                  
##   [9] colorspace_2.0-1              highr_0.8                    
##  [11] knitr_1.31                    rstudioapi_0.13              
##  [13] ggsignif_0.6.0                labeling_0.4.2               
##  [15] GenomeInfoDbData_1.2.1        bit64_4.0.5                  
##  [17] farver_2.0.3                  vctrs_0.3.8                  
##  [19] generics_0.1.0                xfun_0.20                    
##  [21] R6_2.5.0                      locfit_1.5-9.4               
##  [23] manipulateWidget_0.10.1       bitops_1.0-6                 
##  [25] cachem_1.0.1                  assertthat_0.2.1             
##  [27] promises_1.1.1                nnet_7.3-15                  
##  [29] gtable_0.3.0                  affy_1.62.0                  
##  [31] klippy_0.0.0.9500             rlang_0.4.10                 
##  [33] splines_3.6.3                 rstatix_0.6.0                
##  [35] broom_0.7.6                   checkmate_2.0.0              
##  [37] BiocManager_1.30.10           yaml_2.2.1                   
##  [39] abind_1.4-5                   modelr_0.1.8                 
##  [41] crosstalk_1.1.1               backports_1.2.1              
##  [43] httpuv_1.5.5                  Hmisc_4.4-2                  
##  [45] tools_3.6.3                   affyio_1.54.0                
##  [47] ellipsis_0.3.2                base64enc_0.1-3              
##  [49] progress_1.2.2                zlibbioc_1.30.0              
##  [51] RCurl_1.98-1.2                prettyunits_1.1.1            
##  [53] rpart_4.1-15                  cowplot_1.1.1                
##  [55] haven_2.3.1                   cluster_2.1.2                
##  [57] fs_1.5.0                      data.table_1.13.6            
##  [59] RSpectra_0.16-0               openxlsx_4.2.3               
##  [61] reprex_2.0.0                  hms_1.0.0                    
##  [63] mime_0.9                      evaluate_0.14                
##  [65] xtable_1.8-4                  XML_3.99-0.3                 
##  [67] rio_0.5.16                    jpeg_0.1-8.1                 
##  [69] readxl_1.3.1                  gridExtra_2.3                
##  [71] compiler_3.6.3                ellipse_0.4.2                
##  [73] crayon_1.4.1                  htmltools_0.5.1.1            
##  [75] later_1.1.0.1                 Formula_1.2-4                
##  [77] geneplotter_1.62.0            lubridate_1.7.10             
##  [79] DBI_1.1.1                     rappdirs_0.3.2               
##  [81] Matrix_1.3-2                  car_3.0-10                   
##  [83] cli_2.5.0                     pkgconfig_2.0.3              
##  [85] foreign_0.8-75                xml2_1.3.2                   
##  [87] rARPACK_0.11-0                annotate_1.62.0              
##  [89] webshot_0.5.2                 XVector_0.24.0               
##  [91] rvest_1.0.0                   digest_0.6.27                
##  [93] rmarkdown_2.6                 cellranger_1.1.0             
##  [95] fastmatch_1.1-0               htmlTable_2.1.0              
##  [97] shiny_1.6.0                   lifecycle_1.0.0              
##  [99] jsonlite_1.7.2                carData_3.0-4                
## [101] limma_3.40.6                  fansi_0.4.2                  
## [103] pillar_1.6.1                  lattice_0.20-41              
## [105] fastmap_1.1.0                 httr_1.4.2                   
## [107] survival_3.2-7                interactiveDisplayBase_1.22.0
## [109] glue_1.4.2                    zip_2.1.1                    
## [111] png_0.1-7                     bit_4.0.4                    
## [113] stringi_1.5.3                 blob_1.2.1                   
## [115] latticeExtra_0.6-29           memoise_2.0.0

 

---

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu